Roey Angel 2021-05-21
This analysis explores the taxonomical ditribution patters in the different samples, based on the DADA2-produced sequences. Large parts of this script are based on this protocol and the accompanying publichation by Callahan and colleagues (2016).
set.seed(1000)
min_lib_size <- 5000
metadata_path <- "./"
data_path <- "./DADA2_pseudo/"
Metadata_table <- "TeaTime_joint_Fungi_metadata.csv"
Seq_table <- "DADA2.seqtab_nochim.tsv"
Tax_table <- "DADA2.taxa_unite.tsv"
Proj_name <- "TeaTime4Schools"Read abundance table, taxonomic classification and metadata into a phyloseq object. Also remove sequences detected as contaminants in 04_Decontamination.html.
# read OTU mat from data file
Raw_data <- read_tsv(paste0(data_path, Seq_table),
trim_ws = TRUE)
contaminant_seqs <- read_csv(paste0(data_path, "decontam_contaminants.csv"),
trim_ws = TRUE,
col_names = FALSE)
Raw_data %<>% # remove contaminant OTUs.
# .[, -grep("CTRL", colnames(.))] %>% # remove ext. cont.
.[!(Raw_data$`#OTU` %in% contaminant_seqs$X1), ]
Raw_data[, 2:ncol(Raw_data)] %>%
t() %>%
as.data.frame() -> abundance_mat # convert to abundance matrix
colnames(abundance_mat) <- pull(Raw_data, "#OTU") # add sequence names
# Read metadata file
read_csv(paste0(metadata_path, Metadata_table),
trim_ws = TRUE) %>%
mutate_at(
c(
"Workshop",
"Season",
"Run",
"Type",
"Sample type",
"Field",
"Replicate",
"Control",
"Gene"
),
~factor(.)
) %>%
mutate_at(c("Extr. Date", "PCR products_ITS_send for seq"), ~as.Date(., "%d.%m.%Y")) ->
Metadata
Metadata$Season %<>% fct_relevel("Winter", "Spring", "Summer", "Autumn")
Metadata$Read1_file <- str_replace(Metadata$Read1_file, "(.*)_L001_R1_001.fastq.gz|\\.1\\.fastq.gz", "\\1")
Metadata <- Metadata[Metadata$Read1_file %in% rownames(abundance_mat), ] # remove metadata rows if the samples did not go through qual processing
# Order abundance_mat samples according to the metadata
sample_order <- match(rownames(abundance_mat), Metadata$Read1_file)
# Order abundance_mat samples according to the metadata
sample_order <- match(rownames(abundance_mat), Metadata$Read1_file)
abundance_mat %<>%
rownames_to_column('sample_name') %>%
arrange(., sample_order) %>%
column_to_rownames('sample_name') # needed for phyloseq
Metadata$Library.size <- rowSums(abundance_mat)
Metadata <- data.frame(row.names = Metadata$Read1_file, Metadata)
# read taxonomy from data file
Raw_tax_data <- read_tsv(paste0(data_path, Tax_table),
trim_ws = TRUE, col_names = TRUE)
Raw_tax_data %<>%
.[!(Raw_tax_data$`Seq#` %in% contaminant_seqs$X1),] %>% # remove contaminant OTUs.
mutate_all(funs(replace(., is.na(.), "Unclassified"))) # I think mutaute_all is unnecessary here because replace(., is.na(.), "Unclassified") alone should work
Raw_tax_data %>%
select(.,
Kingdom = V8,
Phylum = V9,
Class = V10,
Order = V11,
Family = V12,
Genus = V13) %>%
cbind(Name = colnames(abundance_mat),. ) ->
Taxonomy.bs
Raw_tax_data %>%
select(.,
Kingdom,
Phylum,
Class,
Order,
Family,
Genus) %>%
map_dfr(~str_remove_all(., "^.__")) %>% # remove "k__" prefixes
add_column(Seq = colnames(abundance_mat)) %>%
column_to_rownames("Seq") %>%
set_colnames(c("Domain", "Phylum", "Class", "Order", "Family", "Genus")) %>%
as.matrix() -> # because tax_table() needs a matrix for retaining row.names
Taxonomy
# generate phyloseq object
Ps_obj <- phyloseq(otu_table(abundance_mat, taxa_are_rows = FALSE),
tax_table(Taxonomy),
sample_data(Metadata)
)
Ps_obj <- filter_taxa(Ps_obj, function(x) sum(x) > 0, TRUE) # remove 0 abundance taxa
Ps_obj <- subset_samples(Ps_obj, sample_sums(Ps_obj) > 0) # remove 0 abundance samples
# Remove mock and control samples
Ps_obj <- subset_samples(Ps_obj, Type != "Control" & Type != "Mock")
# Create a grouping variable for merging
sample_data(Ps_obj) %<>%
as(., "data.frame") %>%
# get_variable(., c("Sample.type", "Field", "Season", "Replicate")) %>%
unite(., "Description", c("Sample.type", "Field", "Season", "Replicate"), remove = FALSE)
sample_data(Ps_obj)$Description %<>% as_factor(.)
# merged_Ps_obj <- merge_samples(Ps_obj, "Description")
# merged_SD <- merge_samples(sample_data(Ps_obj), "Description")
Ps_obj_merged <- MergeSamples(Ps_obj, grouping_name = "Description")First, let’s look at the taxonomic distribution
table(tax_table(Ps_obj_merged)[, "Domain"], exclude = NULL)##
## Fungi
## 3620
table(tax_table(Ps_obj_merged)[, "Class"], exclude = NULL)##
## Agaricomycetes Agaricostilbomycetes
## 289 17
## Archaeosporomycetes Basidiobolomycetes
## 2 1
## Blastocladiomycetes Classiculomycetes
## 3 1
## Cystobasidiomycetes Dothideomycetes
## 10 336
## Eurotiomycetes Geminibasidiomycetes
## 120 5
## Glomeromycetes Kickxellomycetes
## 97 13
## Laboulbeniomycetes Lecanoromycetes
## 9 17
## Leotiomycetes Malasseziomycetes
## 153 1
## Microbotryomycetes Mortierellomycetes
## 62 67
## Mucoromycetes Olpidiomycetes
## 6 17
## Orbiliomycetes Paraglomeromycetes
## 29 1
## Pezizomycetes Pucciniomycetes
## 95 2
## Rhizophlyctidomycetes Rhizophydiomycetes
## 9 15
## Rozellomycotina_cls_Incertae_sedis Saccharomycetes
## 7 25
## Sordariomycetes Spizellomycetes
## 607 17
## Taphrinomycetes Tremellomycetes
## 6 122
## Unclassified Ustilaginomycetes
## 1441 1
## Wallemiomycetes Zoopagomycetes
## 1 16
# table(tax_table(Ps_obj_merged)[, "Family"], exclude = NULL)Now let’s remove some taxa which are obvious artefacts (e.g. unclassified domain)
domains2remove <- c("", "Bacteria", "Archaea", "Unclassified")
orders2remove <- c("")
families2remove <- c("")
Ps_obj_filt <- subset_taxa(Ps_obj_merged, !is.na(Phylum) &
!Domain %in% domains2remove &
!Order %in% orders2remove &
!Family %in% families2remove)
Taxonomy.bs %<>%
filter(Taxonomy.bs$Name %in% row.names(Ps_obj_filt@tax_table))
sample_data(Ps_obj_filt)$Lib.size <- rowSums(otu_table(Ps_obj_filt))Now let’s explore the prevalence of different taxa in the database. Prevalence is the number of samples in which a taxa appears at least once. So “Mean prevalence” refers to in how many samples does a sequence belonging to the phylum appears on average, and “Sum prevalence” is the sum of all samples where any sequence from the phylum appears,
prevdf <- apply(X = otu_table(Ps_obj_filt),
MARGIN = ifelse(taxa_are_rows(Ps_obj_filt), yes = 1, no = 2),
FUN = function(x){sum(x > 0)})
# Add taxonomy and total read counts to this data.frame
prevdf <- data.frame(Prevalence = prevdf,
TotalAbundance = taxa_sums(Ps_obj_filt),
tax_table(Ps_obj_filt))
prevdf %>%
group_by(Phylum) %>%
summarise(`Mean prevalence` = mean(Prevalence),
`Sum prevalence` = sum(Prevalence)) ->
Prevalence_phylum_summary
Prevalence_phylum_summary %>%
kable(., digits = c(0, 1, 0)) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = F)| Phylum | Mean prevalence | Sum prevalence |
|---|---|---|
| Ascomycota | 6.5 | 13291 |
| Basidiobolomycota | 4.0 | 4 |
| Basidiomycota | 4.7 | 3050 |
| Blastocladiomycota | 3.7 | 11 |
| Chytridiomycota | 3.2 | 246 |
| Glomeromycota | 3.1 | 324 |
| Kickxellomycota | 1.2 | 16 |
| Mortierellomycota | 9.8 | 775 |
| Mucoromycota | 2.3 | 14 |
| Olpidiomycota | 8.2 | 140 |
| Rozellomycota | 4.2 | 101 |
| Unclassified | 2.6 | 1472 |
| Zoopagomycota | 4.0 | 89 |
prevdf %>%
group_by(Order) %>%
summarise(`Mean prevalence` = mean(Prevalence),
`Sum prevalence` = sum(Prevalence)) ->
Prevalence_order_summary
Prevalence_order_summary %>%
kable(., digits = c(0, 1, 0)) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = F)| Order | Mean prevalence | Sum prevalence |
|---|---|---|
| Agaricales | 3.1 | 489 |
| Agaricostilbales | 5.7 | 97 |
| Archaeosporales | 1.0 | 2 |
| Auriculariales | 5.1 | 71 |
| Basidiobolales | 4.0 | 4 |
| Blastocladiales | 3.7 | 11 |
| Botryosphaeriales | 14.5 | 29 |
| Branch06 | 3.2 | 26 |
| Cantharellales | 3.2 | 131 |
| Capnodiales | 15.6 | 219 |
| Chaetosphaeriales | 2.0 | 6 |
| Chaetothyriales | 5.1 | 318 |
| Classiculales | 1.0 | 1 |
| Coniochaetales | 1.7 | 5 |
| Corticiales | 2.8 | 25 |
| Cystobasidiales | 1.0 | 3 |
| Cystobasidiomycetes\_ord\_Incertae\_sedis | 4.3 | 13 |
| Cystofilobasidiales | 18.6 | 522 |
| Diaporthales | 1.0 | 3 |
| Diversisporales | 5.0 | 15 |
| Dothideales | 13.0 | 26 |
| Dothideomycetes\_ord\_Incertae\_sedis | 7.0 | 28 |
| Erysiphales | 4.0 | 12 |
| Erythrobasidiales | 1.0 | 2 |
| Eurotiales | 6.4 | 264 |
| Filobasidiales | 7.6 | 198 |
| Geastrales | 1.0 | 1 |
| Geminibasidiales | 2.0 | 10 |
| Glomerales | 3.0 | 285 |
| Glomerellales | 6.6 | 146 |
| GS02 | 2.5 | 5 |
| GS05 | 3.0 | 3 |
| GS07 | 6.8 | 27 |
| Helotiales | 12.3 | 1667 |
| Holtermanniales | 13.5 | 81 |
| Hymenochaetales | 1.7 | 5 |
| Hypocreales | 12.7 | 2845 |
| Kickxellales | 1.2 | 16 |
| Leucosporidiales | 12.3 | 111 |
| Magnaporthales | 5.0 | 35 |
| Malasseziales | 1.0 | 1 |
| Melanosporales | 2.8 | 14 |
| Microascales | 3.6 | 64 |
| Microbotryomycetes\_ord\_Incertae\_sedis | 4.9 | 132 |
| Mortierellales | 11.1 | 744 |
| Mucorales | 2.3 | 14 |
| Myrmecridiales | 1.0 | 2 |
| Olpidiales | 8.2 | 140 |
| Onygenales | 2.1 | 35 |
| Orbiliales | 4.4 | 111 |
| Paraglomerales | 8.0 | 8 |
| Pezizales | 4.6 | 433 |
| Phacidiales | 9.0 | 9 |
| Phomatosporales | 7.0 | 7 |
| Platygloeales | 1.5 | 3 |
| Pleosporales | 7.1 | 2123 |
| Pleurotheciales | 4.4 | 31 |
| Polyporales | 2.0 | 20 |
| Pyxidiophorales | 11.6 | 104 |
| Rhizophlyctidales | 3.6 | 32 |
| Rhizophydiales | 3.3 | 49 |
| Russulales | 7.0 | 98 |
| Saccharomycetales | 3.6 | 89 |
| Sebacinales | 5.0 | 114 |
| Sordariales | 7.1 | 1285 |
| Spizellomycetales | 1.6 | 27 |
| Sporidiobolales | 5.0 | 109 |
| Taphrinales | 5.0 | 30 |
| Thelebolales | 3.4 | 31 |
| Thelephorales | 1.0 | 1 |
| Togniniales | 1.0 | 3 |
| Trechisporales | 2.3 | 7 |
| Tremellales | 7.5 | 414 |
| Trichosphaeriales | 1.0 | 2 |
| Trichosporonales | 8.8 | 35 |
| Tubeufiales | 9.3 | 28 |
| Unclassified | 3.1 | 4869 |
| Ustilaginales | 1.0 | 1 |
| Venturiales | 4.7 | 28 |
| Wallemiales | 1.0 | 1 |
| Xylariales | 8.1 | 459 |
| Zoopagales | 4.3 | 69 |
Based on that I’ll remove all phyla with a prevalence of under 10
Prevalence_phylum_summary %>%
filter(`Sum prevalence` < 10) %>%
select(Phylum) %>%
map(as.character) %>%
unlist() ->
filterPhyla
Ps_obj_filt2 <- subset_taxa(Ps_obj_filt, !Phylum %in% filterPhyla)
Taxonomy.bs %<>%
filter(Taxonomy.bs$Name %in% row.names(Ps_obj_filt2@tax_table))
sample_data(Ps_obj_filt2)$Lib.size <- rowSums(otu_table(Ps_obj_filt2))
print(Ps_obj_filt)## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 3620 taxa and 123 samples ]
## sample_data() Sample Data: [ 123 samples by 26 sample variables ]
## tax_table() Taxonomy Table: [ 3620 taxa by 6 taxonomic ranks ]
print(Ps_obj_filt2)## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 3619 taxa and 123 samples ]
## sample_data() Sample Data: [ 123 samples by 26 sample variables ]
## tax_table() Taxonomy Table: [ 3619 taxa by 6 taxonomic ranks ]
Plot general prevalence features of the phyla
# Subset to the remaining phyla
prevdf_phylum_filt <- subset(prevdf, Phylum %in% get_taxa_unique(Ps_obj_filt2, "Phylum"))
ggplot(prevdf_phylum_filt,
aes(TotalAbundance, Prevalence / nsamples(Ps_obj_filt2), color = Phylum)) +
# Include a guess for parameter
geom_hline(yintercept = 0.05,
alpha = 0.5,
linetype = 2) + geom_point(size = 2, alpha = 0.7) +
scale_x_log10() + xlab("Total Abundance") + ylab("Prevalence [Frac. Samples]") +
facet_wrap( ~ Phylum) + theme(legend.position = "none")Plot general prevalence features of the top 20 orders
# Subset to the remaining phyla
prevdf_order_filt <- subset(prevdf, Order %in% get_taxa_unique(Ps_obj_filt2, "Order"))
# grab the top 30 most abundant orders
prevdf_order_filt %>%
group_by(Order) %>%
summarise(Combined.abundance = sum(TotalAbundance)) %>%
arrange(desc(Combined.abundance)) %>%
.[1:30, "Order"] ->
Orders2plot
prevdf_order_filt2 <- subset(prevdf, Order %in% Orders2plot$Order)
ggplot(prevdf_order_filt2,
aes(TotalAbundance, Prevalence / nsamples(Ps_obj_filt2), color = Order)) +
# Include a guess for parameter
geom_hline(yintercept = 0.05,
alpha = 0.5,
linetype = 2) + geom_point(size = 2, alpha = 0.7) +
scale_x_log10() + xlab("Total Abundance") + ylab("Prevalence [Frac. Samples]") +
facet_wrap( ~ Order) + theme(legend.position = "none")Plot general prevalence features of the phyla excluding soil samples
Ps_obj_filt_Teas <- subset_samples(Ps_obj_filt2, Type != "Soil")
prevdf_teas <- apply(X = otu_table(Ps_obj_filt_Teas),
MARGIN = ifelse(taxa_are_rows(Ps_obj_filt_Teas), yes = 1, no = 2),
FUN = function(x){sum(x > 0)})
# Add taxonomy and total read counts to this data.frame
prevdf_teas <- data.frame(Prevalence = prevdf_teas,
TotalAbundance = taxa_sums(Ps_obj_filt_Teas),
tax_table(Ps_obj_filt_Teas))
# Subset to the remaining phyla
prevdf_phylum_filt <- subset(prevdf_teas, Phylum %in% get_taxa_unique(Ps_obj_filt_Teas, "Phylum"))
ggplot(prevdf_phylum_filt,
aes(TotalAbundance, Prevalence / nsamples(Ps_obj_filt_Teas), color = Phylum)) +
# Include a guess for parameter
geom_hline(yintercept = 0.05,
alpha = 0.5,
linetype = 2) + geom_point(size = 2, alpha = 0.7) +
scale_x_log10() + xlab("Total Abundance") + ylab("Prevalence [Frac. Samples]") +
facet_wrap( ~ Phylum) + theme(legend.position = "none")Plot general prevalence features of the top 20 orders (teabag samples only)
# Subset to the remaining phyla
prevdf_order_filt <- subset(prevdf_teas, Order %in% get_taxa_unique(Ps_obj_filt_Teas, "Order"))
# grab the top 30 most abundant orders
prevdf_order_filt %>%
group_by(Order) %>%
summarise(Combined.abundance = sum(TotalAbundance)) %>%
arrange(desc(Combined.abundance)) %>%
.[1:30, "Order"] ->
Orders2plot
prevdf_order_filt2 <- subset(prevdf, Order %in% Orders2plot$Order)
ggplot(prevdf_order_filt2,
aes(TotalAbundance, Prevalence / nsamples(Ps_obj_filt_Teas), color = Order)) +
# Include a guess for parameter
geom_hline(yintercept = 0.05,
alpha = 0.5,
linetype = 2) + geom_point(size = 2, alpha = 0.7) +
scale_x_log10() + xlab("Total Abundance") + ylab("Prevalence [Frac. Samples]") +
facet_wrap( ~ Order) + theme(legend.position = "none") #### Unsupervised filtering by prevalence I’ll remove all sequences which appear in less than 5% of the samples
# Define prevalence threshold as 5% of total samples
prevalenceThreshold <- 0.05 * nsamples(Ps_obj_filt)
prevalenceThreshold## [1] 6.15
# Execute prevalence filter, using `prune_taxa()` function
keepTaxa <-
row.names(prevdf_phylum_filt)[(prevdf_phylum_filt$Prevalence >= prevalenceThreshold)]
Ps_obj_filt3 <- prune_taxa(keepTaxa, Ps_obj_filt2)
Taxonomy.bs %<>%
filter(Taxonomy.bs$Name %in% row.names(Ps_obj_filt3@tax_table))
sample_data(Ps_obj_filt3)$Lib.size <- rowSums(otu_table(Ps_obj_filt3))
print(Ps_obj_filt2)## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 3619 taxa and 123 samples ]
## sample_data() Sample Data: [ 123 samples by 26 sample variables ]
## tax_table() Taxonomy Table: [ 3619 taxa by 6 taxonomic ranks ]
print(Ps_obj_filt3)## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 402 taxa and 123 samples ]
## sample_data() Sample Data: [ 123 samples by 26 sample variables ]
## tax_table() Taxonomy Table: [ 402 taxa by 6 taxonomic ranks ]
This removed 3217 or 89% of the ESVs!.
However all these removed ESVs accounted for only:
prevdf_phylum_filt %>%
arrange(., Prevalence) %>%
group_by(Prevalence > prevalenceThreshold) %>%
summarise(Abundance = sum(TotalAbundance)) %>%
mutate(`Rel. Ab.` = percent(Abundance / sum(Abundance))) %>%
kable(., digits = c(0, 1, 0)) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = F)| Prevalence > prevalenceThreshold | Abundance | Rel. Ab. |
|---|---|---|
| FALSE | 290789 | 11% |
| TRUE | 2315738 | 89% |
So it’s fine to remove them.
Let’s make a fasta file after filtering the sequences:
Seq_file <- "DADA2.Seqs.fa"
readDNAStringSet(paste0(data_path, Seq_file)) %>%
.[taxa_names(Ps_obj_filt3)] %>%
writeXStringSet(., filepath = paste0(data_path, str_remove(Seq_file, ".fa*"), "_filtered.fa"), format = "fasta", width = 1000)Plot general prevalence features of the phyla after filtering
# Subset to the remaining phyla
prevdf_phylum_filt <- subset(prevdf, Phylum %in% get_taxa_unique(Ps_obj_filt3, "Phylum"))
ggplot(prevdf_phylum_filt,
aes(TotalAbundance, Prevalence / nsamples(Ps_obj_filt3), color = Phylum)) +
# Include a guess for parameter
geom_hline(yintercept = 0.05,
alpha = 0.5,
linetype = 2) + geom_point(size = 2, alpha = 0.7) +
scale_x_log10() + xlab("Total Abundance") + ylab("Prevalence [Frac. Samples]") +
facet_wrap( ~ Phylum) + theme(legend.position = "none")Plot general prevalence features of the top 20 orders
# Subset to the remaining phyla
prevdf_order_filt <- subset(prevdf, Order %in% get_taxa_unique(Ps_obj_filt3, "Order"))
# grab the top 30 most abundant orders
prevdf_order_filt %>%
group_by(Order) %>%
summarise(Combined.abundance = sum(TotalAbundance)) %>%
arrange(desc(Combined.abundance)) %>%
.[1:30, "Order"] ->
Orders2plot
prevdf_order_filt2 <- subset(prevdf, Order %in% Orders2plot$Order)
ggplot(prevdf_order_filt2,
aes(TotalAbundance, Prevalence / nsamples(Ps_obj_filt3), color = Order)) +
# Include a guess for parameter
geom_hline(yintercept = 0.05,
alpha = 0.5,
linetype = 2) + geom_point(size = 2, alpha = 0.7) +
scale_x_log10() + xlab("Total Abundance") + ylab("Prevalence [Frac. Samples]") +
facet_wrap( ~ Order) + theme(legend.position = "none")Plot general prevalence features of the phyla excluding soil samples
# Subset to the remaining phyla
prevdf_phylum_filt_teas <- subset(prevdf_teas, Phylum %in% get_taxa_unique(Ps_obj_filt_Teas, "Phylum"))
ggplot(prevdf_phylum_filt,
aes(TotalAbundance, Prevalence / nsamples(Ps_obj_filt_Teas), color = Phylum)) +
# Include a guess for parameter
geom_hline(yintercept = 0.05,
alpha = 0.5,
linetype = 2) + geom_point(size = 2, alpha = 0.7) +
scale_x_log10() + xlab("Total Abundance") + ylab("Prevalence [Frac. Samples]") +
facet_wrap( ~ Phylum) + theme(legend.position = "none")Plot general prevalence features of the top 20 orders
# Subset to the remaining phyla
prevdf_order_filt <- subset(prevdf_teas, Order %in% get_taxa_unique(Ps_obj_filt_Teas, "Order"))
# grab the top 30 most abundant orders
prevdf_order_filt %>%
group_by(Order) %>%
summarise(Combined.abundance = sum(TotalAbundance)) %>%
arrange(desc(Combined.abundance)) %>%
.[1:30, "Order"] ->
Orders2plot
prevdf_order_filt2 <- subset(prevdf, Order %in% Orders2plot$Order)
ggplot(prevdf_order_filt2,
aes(TotalAbundance, Prevalence / nsamples(Ps_obj_filt_Teas), color = Order)) +
# Include a guess for parameter
geom_hline(yintercept = 0.05,
alpha = 0.5,
linetype = 2) + geom_point(size = 2, alpha = 0.7) +
scale_x_log10() + xlab("Total Abundance") + ylab("Prevalence [Frac. Samples]") +
facet_wrap( ~ Order) + theme(legend.position = "none")# remove samples with reads < min_lib_size
Ps_obj_filt3 %>%
subset_samples(., sample_sums(Ps_obj_filt3) > min_lib_size) %>%
filter_taxa(., function(x)
sum(x) > 0, TRUE) ->
# transform_sample_counts(., function(x)
# x / sum(x) * 100) ->
Ps_obj_filt3_subset
Ps_obj_filt3_subset %>%
get_variable() %>%
mutate_at(., "Sample.type",
~fct_relevel(., levels = c("Soil", "Green tea", "Rooibos"))) %>%
mutate_at(., "Season",
~fct_relevel(., levels = c("Winter", "Spring", "Summer", "Autumn"))) %>%
arrange(Field, Sample.type, Season, Replicate) %>%
pull(Description) %>%
as.character() ->
Sample.order
Ps_obj_filt3_subset_ra <- transform_sample_counts(Ps_obj_filt3_subset, function(x){x / sum(x)} * 100)
# grabthe top 100 most abundant OTUs
Ps_obj_filt3_subset_100 <-
prune_taxa(names(sort(taxa_sums(Ps_obj_filt3_subset_ra), TRUE)[1:100]), Ps_obj_filt3_subset_ra)
plot_heatmap(
Ps_obj_filt3_subset_100,
method = NULL,
distance = NULL,
sample.label = "Description",
taxa.label = "Order",
sample.order = Sample.order,
low = "#000033",
high = "#FF3300"
)Ps_obj_filt3_subset_ra %>%
subset_taxa(., Phylum == "Ascomycota") %>%
transform_sample_counts(., function(x) x / sum(x) * 100) ->
Ps_obj_filt3_subset_asco
plot_heatmap(
Ps_obj_filt3_subset_asco,
method = NULL,
distance = NULL,
sample.label = "Description",
sample.order = Sample.order,
taxa.label = "Order",
low = "#000033",
high = "#FF3300"
)Ps_obj_filt3_subset_ra %>%
subset_taxa(., Phylum == "Basidiomycota") %>%
transform_sample_counts(., function(x) x / sum(x) * 100) ->
Ps_obj_filt3_subset_basido
plot_heatmap(
Ps_obj_filt3_subset_basido,
method = NULL,
distance = NULL,
sample.label = "Description",
sample.order = Sample.order,
taxa.label = "Order",
low = "#000033",
high = "#FF3300"
)Ps_obj_filt3_subset_ra %>%
subset_taxa(., Phylum == "Chytridiomycota") %>%
transform_sample_counts(., function(x) x / sum(x) * 100) ->
Ps_obj_filt3_subset_chytri
plot_heatmap(
Ps_obj_filt3_subset_chytri,
method = NULL,
distance = NULL,
sample.label = "Description",
sample.order = Sample.order,
taxa.label = "Order",
low = "#000033",
high = "#FF3300"
) Ps_obj_filt3_subset_ra %>%
subset_taxa(., Phylum == "Mortierellomycota") %>%
transform_sample_counts(., function(x) x / sum(x) * 100) ->
Ps_obj_filt3_subset_morti
plot_heatmap(
Ps_obj_filt3_subset_morti,
method = NULL,
distance = NULL,
sample.label = "Description",
sample.order = Sample.order,
taxa.label = "Phylum",
low = "#000033",
high = "#FF3300"
)Ps_obj_filt3_subset_ra %>%
subset_taxa(., Order == "Helotiales") %>%
transform_sample_counts(., function(x) x / sum(x) * 100) ->
Ps_obj_filt3_subset_helo
plot_heatmap(
Ps_obj_filt3_subset_helo,
method = NULL,
distance = NULL,
sample.label = "Description",
sample.order = Sample.order,
taxa.label = "Order",
low = "#000033",
high = "#FF3300"
)Ps_obj_filt3_subset_ra %>%
subset_taxa(., Order == "Hypocreales") %>%
transform_sample_counts(., function(x) x / sum(x) * 100) ->
Ps_obj_filt3_subset_hypo
plot_heatmap(
Ps_obj_filt3_subset_hypo,
method = NULL,
distance = NULL,
sample.label = "Description",
sample.order = Sample.order,
taxa.label = "Order",
low = "#000033",
high = "#FF3300"
)Ps_obj_filt3_subset_ra %>%
subset_taxa(., Order == "Pleosporales") %>%
transform_sample_counts(., function(x) x / sum(x) * 100) ->
Ps_obj_filt3_subset_pleo
plot_heatmap(
Ps_obj_filt3_subset_pleo,
method = NULL,
distance = NULL,
sample.label = "Description",
sample.order = Sample.order,
taxa.label = "Order",
low = "#000033",
high = "#FF3300"
)Let’s look at the agglomerated taxa
Ps_obj_filt3_subset_ra %>%
tax_glom(., "Order", NArm = TRUE) ->
Ps_obj_filt_glom
plot_heatmap(
Ps_obj_filt_glom,
# method = "NMDS",
# distance = "bray",
sample.order = Sample.order,
sample.label = "Description",
taxa.label = "Order",
taxa.order = "Order",
low = "#000033",
high = "#FF3300"
) + theme_bw(base_size = 20) + theme(axis.text.x = element_text(hjust = 1.0, angle = 90.0))PlotAbundance(Ps_obj_filt3_subset_ra, phylum2plot = "Ascomycota")# plotBefore <- PlotAbundance(Ps_obj_filt3_subset, taxa = "Proteobacteria")
# plotAfter <- PlotAbundance(Ps_obj_filt3_subset_ra, taxa = "Proteobacteria")
# Combine each plot into one graphic.
# grid.arrange(nrow = 2, plotBefore, plotAfter)Ps_obj_filt3_subset_ra_taxa <- subset_taxa(Ps_obj_filt3_subset_ra, Order == "Helotiales")
PlotAbundance(Ps_obj_filt3_subset_ra_taxa, phylum2plot = "Ascomycota", Facet = "Genus")Ps_obj_filt3_subset_df <- psmelt(Ps_obj_filt3_subset)
# Ps_obj_filt3_subset_df %>%
# filter(Species == "Epibolus pulchripes") ->
# Ps_obj_filt3_subset_df_EP
# PlotViolin(Ps_obj_filt3_subset_df_EP)
PlotViolin(Ps_obj_filt3_subset_df)Ps_obj_filt3_glom <- tax_glom(Ps_obj_filt3_subset_ra,
"Order",
NArm = TRUE)
Ps_obj_filt3_glom_DF <- psmelt(Ps_obj_filt3_glom)
Ps_obj_filt3_glom_DF$Order %<>% as.character()
# group dataframe by Phylum, calculate median rel. abundance
Ps_obj_filt3_glom_DF %>%
group_by(Order) %>%
summarise(median = median(Abundance)) ->
medians
# find Phyla whose rel. abund. is less than 0.5%
Rare_phyla <- medians[medians$median <= 0.005, ]$Order
# change their name to "Rare"
Ps_obj_filt3_glom_DF[Ps_obj_filt3_glom_DF$Order %in% Rare_phyla, ]$Order <- 'Rare'
# re-group
Ps_obj_filt3_glom_DF %>%
group_by(Sample, Sample.name, Order, Sample.type, Field) %>%
summarise(Abundance = sum(Abundance)) ->
Ps_obj_filt3_glom_DF_2plot
# ab.taxonomy$Freq <- sqrt(ab.taxonomy$Freq)
# Ps_obj_filt3_glom_rel_DF$Phylum %<>% sub("unclassified", "Unclassified", .)
# Ps_obj_filt3_glom_rel_DF$Phylum %<>% sub("uncultured", "Unclassified", .)
Ps_obj_filt3_glom_DF_2plot %>%
group_by(Sample) %>%
filter(Order == "Rare") %>%
summarise(`Rares (%)` = sum(Abundance * 100)) ->
Rares
# Percentage of reads classified as rare
Rares %>%
kable(., digits = c(2), caption = "Percentage of reads per sample classified as rare:") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = F)| Sample | Rares (%) |
|---|---|
| Green\_tea\_Auboden\_Autumn\_1 | 2556.79 |
| Green\_tea\_Auboden\_Autumn\_2 | 1478.05 |
| Green\_tea\_Auboden\_Autumn\_3 | 155.19 |
| Green\_tea\_Auboden\_Autumn\_4 | 475.83 |
| Green\_tea\_Auboden\_Spring\_1 | 47.91 |
| Green\_tea\_Auboden\_Spring\_2 | 123.91 |
| Green\_tea\_Auboden\_Spring\_3 | 45.19 |
| Green\_tea\_Auboden\_Spring\_4 | 13.05 |
| Green\_tea\_Auboden\_Summer\_1 | 101.25 |
| Green\_tea\_Auboden\_Summer\_2 | 17.84 |
| Green\_tea\_Auboden\_Summer\_3 | 229.51 |
| Green\_tea\_Auboden\_Summer\_4 | 61.60 |
| Green\_tea\_Auboden\_Winter\_1 | 190.33 |
| Green\_tea\_Auboden\_Winter\_2 | 632.67 |
| Green\_tea\_Auboden\_Winter\_3 | 378.26 |
| Green\_tea\_Auboden\_Winter\_4 | 502.86 |
| Green\_tea\_Braunerde\_Autumn\_1 | 115.76 |
| Green\_tea\_Braunerde\_Autumn\_2 | 7.59 |
| Green\_tea\_Braunerde\_Autumn\_3 | 26.60 |
| Green\_tea\_Braunerde\_Autumn\_4 | 42.17 |
| Green\_tea\_Braunerde\_Spring\_1 | 138.24 |
| Green\_tea\_Braunerde\_Spring\_2 | 199.20 |
| Green\_tea\_Braunerde\_Spring\_3 | 267.15 |
| Green\_tea\_Braunerde\_Spring\_4 | 251.97 |
| Green\_tea\_Braunerde\_Summer\_1 | 635.72 |
| Green\_tea\_Braunerde\_Summer\_2 | 434.43 |
| Green\_tea\_Braunerde\_Summer\_3 | 313.91 |
| Green\_tea\_Braunerde\_Summer\_4 | 693.79 |
| Green\_tea\_Braunerde\_Winter\_1 | 1428.98 |
| Green\_tea\_Braunerde\_Winter\_2 | 827.83 |
| Green\_tea\_Braunerde\_Winter\_3 | 779.16 |
| Green\_tea\_Braunerde\_Winter\_4 | 1071.17 |
| Green\_tea\_Kolluvisol\_Autumn\_1 | 130.20 |
| Green\_tea\_Kolluvisol\_Autumn\_2 | 28.55 |
| Green\_tea\_Kolluvisol\_Autumn\_3 | 25.35 |
| Green\_tea\_Kolluvisol\_Autumn\_4 | 273.94 |
| Green\_tea\_Kolluvisol\_Spring\_1 | 156.45 |
| Green\_tea\_Kolluvisol\_Spring\_3 | 16.46 |
| Green\_tea\_Kolluvisol\_Spring\_4 | 75.85 |
| Green\_tea\_Kolluvisol\_Summer\_2 | 114.57 |
| Green\_tea\_Kolluvisol\_Summer\_3 | 68.23 |
| Green\_tea\_Kolluvisol\_Winter\_1 | 2052.51 |
| Green\_tea\_Kolluvisol\_Winter\_2 | 946.75 |
| Green\_tea\_Kolluvisol\_Winter\_3 | 1234.56 |
| Green\_tea\_Kolluvisol\_Winter\_4 | 149.13 |
| Green\_tea\_Unburied\_Winter\_1 | 1179.31 |
| Rooibos\_Auboden\_Autumn\_1 | 1285.15 |
| Rooibos\_Auboden\_Autumn\_2 | 1891.21 |
| Rooibos\_Auboden\_Autumn\_3 | 806.06 |
| Rooibos\_Auboden\_Autumn\_4 | 79.47 |
| Rooibos\_Auboden\_Spring\_1 | 97.25 |
| Rooibos\_Auboden\_Spring\_2 | 1720.92 |
| Rooibos\_Auboden\_Spring\_3 | 167.17 |
| Rooibos\_Auboden\_Spring\_4 | 2017.53 |
| Rooibos\_Auboden\_Summer\_1 | 500.39 |
| Rooibos\_Auboden\_Summer\_2 | 310.81 |
| Rooibos\_Auboden\_Summer\_3 | 980.45 |
| Rooibos\_Auboden\_Summer\_4 | 1973.63 |
| Rooibos\_Auboden\_Winter\_1 | 613.94 |
| Rooibos\_Auboden\_Winter\_2 | 297.24 |
| Rooibos\_Auboden\_Winter\_3 | 579.50 |
| Rooibos\_Auboden\_Winter\_4 | 503.64 |
| Rooibos\_Braunerde\_Autumn\_1 | 1244.06 |
| Rooibos\_Braunerde\_Autumn\_2 | 502.69 |
| Rooibos\_Braunerde\_Autumn\_4 | 2624.33 |
| Rooibos\_Braunerde\_Spring\_1 | 878.67 |
| Rooibos\_Braunerde\_Spring\_3 | 4005.45 |
| Rooibos\_Braunerde\_Summer\_1 | 120.75 |
| Rooibos\_Braunerde\_Summer\_2 | 91.24 |
| Rooibos\_Braunerde\_Summer\_4 | 1611.13 |
| Rooibos\_Braunerde\_Winter\_1 | 3547.66 |
| Rooibos\_Braunerde\_Winter\_2 | 1456.62 |
| Rooibos\_Braunerde\_Winter\_3 | 2011.76 |
| Rooibos\_Braunerde\_Winter\_4 | 5062.61 |
| Rooibos\_Kolluvisol\_Autumn\_1 | 1287.79 |
| Rooibos\_Kolluvisol\_Autumn\_2 | 8156.63 |
| Rooibos\_Kolluvisol\_Autumn\_3 | 213.66 |
| Rooibos\_Kolluvisol\_Autumn\_4 | 9357.24 |
| Rooibos\_Kolluvisol\_Spring\_1 | 23.11 |
| Rooibos\_Kolluvisol\_Spring\_3 | 471.42 |
| Rooibos\_Kolluvisol\_Spring\_4 | 1293.25 |
| Rooibos\_Kolluvisol\_Summer\_1 | 196.29 |
| Rooibos\_Kolluvisol\_Summer\_2 | 97.66 |
| Rooibos\_Kolluvisol\_Summer\_4 | 19.07 |
| Rooibos\_Kolluvisol\_Winter\_1 | 3339.25 |
| Rooibos\_Kolluvisol\_Winter\_2 | 2483.24 |
| Rooibos\_Kolluvisol\_Winter\_3 | 1477.80 |
| Rooibos\_Kolluvisol\_Winter\_4 | 782.97 |
| Soil\_Auboden\_Autumn\_1 | 80.41 |
| Soil\_Auboden\_Autumn\_2 | 146.12 |
| Soil\_Auboden\_Autumn\_3 | 63.47 |
| Soil\_Auboden\_Spring\_1 | 268.34 |
| Soil\_Auboden\_Spring\_2 | 131.34 |
| Soil\_Auboden\_Summer\_1 | 179.93 |
| Soil\_Auboden\_Summer\_2 | 73.59 |
| Soil\_Auboden\_Summer\_3 | 131.15 |
| Soil\_Auboden\_Winter\_1 | 551.10 |
| Soil\_Auboden\_Winter\_2 | 200.47 |
| Soil\_Braunerde\_Autumn\_1 | 187.99 |
| Soil\_Braunerde\_Autumn\_2 | 80.04 |
| Soil\_Braunerde\_Autumn\_3 | 403.20 |
| Soil\_Braunerde\_Spring\_1 | 246.67 |
| Soil\_Braunerde\_Spring\_2 | 732.39 |
| Soil\_Braunerde\_Summer\_1 | 438.11 |
| Soil\_Braunerde\_Summer\_2 | 1104.98 |
| Soil\_Braunerde\_Summer\_3 | 408.50 |
| Soil\_Braunerde\_Winter\_1 | 1149.69 |
| Soil\_Braunerde\_Winter\_2 | 322.34 |
| Soil\_Kolluvisol\_Autumn\_1 | 139.00 |
| Soil\_Kolluvisol\_Autumn\_2 | 133.71 |
| Soil\_Kolluvisol\_Autumn\_3 | 469.74 |
| Soil\_Kolluvisol\_Spring\_1 | 187.55 |
| Soil\_Kolluvisol\_Spring\_2 | 82.36 |
| Soil\_Kolluvisol\_Summer\_1 | 174.78 |
| Soil\_Kolluvisol\_Summer\_2 | 136.28 |
| Soil\_Kolluvisol\_Summer\_3 | 216.76 |
| Soil\_Kolluvisol\_Winter\_1 | 2140.22 |
| Soil\_Kolluvisol\_Winter\_2 | 189.77 |
sample_order <- match(Rares$Sample, row.names(get_variable(Ps_obj_filt3_glom)))
Rares %<>% arrange(., sample_order)
Rares %>%
cbind(., get_variable(Ps_obj_filt3_glom)) %>%
group_by(Sample.type) %>%
summarise(`Rares (%)` = mean(`Rares (%)`)) ->
Rares_merged
# Percentage of reads classified as rare
Rares_merged %>%
kable(., digits = c(2), caption = "Percentage of reads per sample classified as rare:") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = F)| Sample.type | Rares (%) |
|---|---|
| Green tea | 449.91 |
| Rooibos | 1575.73 |
| Soil | 359.00 |
Ps_obj_filt3_glom_DF_2plot %>%
group_by(Order) %>%
summarise(sum.Taxa = sum(Abundance)) %>%
arrange(desc(sum.Taxa)) -> Taxa.rank
Ps_obj_filt3_glom_DF_2plot$Order %<>%
factor(., levels = Taxa.rank$Order) %>%
fct_relevel(., "Rare", after = Inf)
p_taxa_box <-
ggplot(Ps_obj_filt3_glom_DF_2plot, aes(x = Order, y = (Abundance))) +
geom_boxplot(aes(group = interaction(Order, Sample.type)), position = position_dodge(width = 0.9), fatten = 1) +
geom_point(
aes(colour = Sample.type),
position = position_jitterdodge(dodge.width = 1),
alpha = 1 / 2,
stroke = 0,
size = 2
) +
# scale_colour_manual(values = pom4, name = "") +
theme_cowplot(font_size = 11, font_family = f_name) +
labs(x = NULL, y = "Relative abundance (%)") +
guides(colour = guide_legend(override.aes = list(size = 5))) +
facet_grid(Field ~ .) +
background_grid(major = "xy",
minor = "none") +
theme(axis.text.x = element_text(
angle = 45,
vjust = 0.9,
hjust = 0.9
))
print(p_taxa_box)dir.create(paste0(data_path, "Krona/"))
for (i in seq(nsamples(Ps_obj_filt3_subset))) {
sample_data <-
data.frame(Abundance = get_taxa(Ps_obj_filt3_subset, i), as(tax_table(Ps_obj_filt3_subset), "matrix"))
sample_data <- sample_data[sample_data[, 1] > 0, ]
write_tsv(sample_data, paste0(data_path, "Krona/", sample_names(Ps_obj_filt3_subset)[i], ".tsv"))
}
sample_data(Ps_obj_filt3_subset)$Sample.type %<>% fct_recode(., Green_tea = "Green tea")
sample_data(Ps_obj_filt3_subset)$Sample.type_field_season <- paste0(
get_variable(Ps_obj_filt3_subset, "Sample.type"),
"_",
get_variable(Ps_obj_filt3_subset, "Field"),
"_",
get_variable(Ps_obj_filt3_subset, "Season")
)
Ps_obj_filt3_subset %>%
phyloseq::merge_samples(., "Sample.type_field_season", fun = mean) ->
merged.Ps_obj
for (i in seq(nsamples(merged.Ps_obj))) {
sample_data <-
data.frame(Abundance = get_taxa(merged.Ps_obj, i), as(tax_table(merged.Ps_obj), "matrix"))
sample_data <- sample_data[sample_data[, 1] > 0, ]
write_tsv(sample_data, paste0(data_path, "Krona/", sample_names(merged.Ps_obj)[i], ".tsv"))
}
list.files(paste0(data_path, "Krona/"), full.names = TRUE) %>%
paste(., collapse = " ") %>%
paste0("/usr/local/bin/ktImportText ", .,
" -o ",
data_path,
"Krona/TeaTime_Krona.html") %>%
system()saveRDS(Ps_obj_filt, file = paste0(data_path, Proj_name, "ITS_filt.RDS"))
saveRDS(Ps_obj_filt3, file = paste0(data_path, Proj_name, "ITS_filt3.RDS"))sessioninfo::session_info() %>%
details::details(
summary = 'Current session info',
open = TRUE
)─ Session info ─────────────────────────────────────────────────────────────────────────
setting value
version R version 4.0.3 (2020-10-10)
os Ubuntu 18.04.5 LTS
system x86_64, linux-gnu
ui X11
language (EN)
collate en_US.UTF-8
ctype en_US.UTF-8
tz Europe/Prague
date 2021-05-21
─ Packages ─────────────────────────────────────────────────────────────────────────────
package * version date lib source
ade4 1.7-16 2020-10-28 [1] CRAN (R 4.0.2)
ape 5.5 2021-04-25 [1] CRAN (R 4.0.3)
assertthat 0.2.1 2019-03-21 [1] CRAN (R 4.0.2)
backports 1.2.1 2020-12-09 [1] CRAN (R 4.0.2)
Biobase 2.50.0 2020-10-27 [1] Bioconductor
BiocGenerics * 0.36.1 2021-04-16 [1] Bioconductor
biomformat 1.18.0 2020-10-27 [1] Bioconductor
Biostrings * 2.58.0 2020-10-27 [1] Bioconductor
broom * 0.7.6 2021-04-05 [1] CRAN (R 4.0.3)
cellranger 1.1.0 2016-07-27 [1] CRAN (R 4.0.2)
cli 2.5.0 2021-04-26 [1] CRAN (R 4.0.3)
clipr 0.7.1 2020-10-08 [1] CRAN (R 4.0.2)
cluster 2.1.2 2021-04-17 [1] CRAN (R 4.0.3)
codetools 0.2-18 2020-11-04 [1] CRAN (R 4.0.2)
colorspace 2.0-1 2021-05-04 [1] CRAN (R 4.0.3)
cowplot * 1.1.1 2020-12-30 [1] CRAN (R 4.0.2)
crayon 1.4.1 2021-02-08 [1] CRAN (R 4.0.3)
data.table 1.14.0 2021-02-21 [1] CRAN (R 4.0.3)
DBI 1.1.1 2021-01-15 [1] CRAN (R 4.0.3)
dbplyr 2.1.1 2021-04-06 [1] CRAN (R 4.0.3)
desc 1.3.0 2021-03-05 [1] CRAN (R 4.0.3)
details 0.2.1 2020-01-12 [1] CRAN (R 4.0.2)
digest 0.6.27 2020-10-24 [1] CRAN (R 4.0.2)
doParallel * 1.0.16 2020-10-16 [1] CRAN (R 4.0.2)
dplyr * 1.0.6 2021-05-05 [1] CRAN (R 4.0.3)
ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.0.3)
evaluate 0.14 2019-05-28 [1] CRAN (R 4.0.2)
extrafont * 0.17 2014-12-08 [1] CRAN (R 4.0.2)
extrafontdb 1.0 2012-06-11 [1] CRAN (R 4.0.2)
fansi 0.4.2 2021-01-15 [1] CRAN (R 4.0.3)
farver 2.1.0 2021-02-28 [1] CRAN (R 4.0.3)
forcats * 0.5.1 2021-01-27 [1] CRAN (R 4.0.3)
foreach * 1.5.1 2020-10-15 [1] CRAN (R 4.0.2)
fs 1.5.0 2020-07-31 [1] CRAN (R 4.0.2)
generics 0.1.0 2020-10-31 [1] CRAN (R 4.0.2)
ggplot2 * 3.3.3 2020-12-30 [1] CRAN (R 4.0.2)
glue 1.4.2 2020-08-27 [1] CRAN (R 4.0.2)
gridExtra * 2.3 2017-09-09 [1] CRAN (R 4.0.2)
gtable 0.3.0 2019-03-25 [1] CRAN (R 4.0.2)
haven 2.4.1 2021-04-23 [1] CRAN (R 4.0.3)
highr 0.9 2021-04-16 [1] CRAN (R 4.0.3)
hms 1.1.0 2021-05-17 [1] CRAN (R 4.0.3)
htmltools 0.5.1.1 2021-01-22 [1] CRAN (R 4.0.3)
httr 1.4.2 2020-07-20 [1] CRAN (R 4.0.2)
igraph 1.2.6 2020-10-06 [1] CRAN (R 4.0.2)
IRanges * 2.24.1 2020-12-12 [1] Bioconductor
iterators * 1.0.13 2020-10-15 [1] CRAN (R 4.0.2)
jsonlite 1.7.2 2020-12-09 [1] CRAN (R 4.0.2)
kableExtra * 1.3.4 2021-02-20 [1] CRAN (R 4.0.3)
knitr * 1.33 2021-04-24 [1] CRAN (R 4.0.3)
labeling 0.4.2 2020-10-20 [1] CRAN (R 4.0.2)
lattice * 0.20-44 2021-05-02 [1] CRAN (R 4.0.3)
lifecycle 1.0.0 2021-02-15 [1] CRAN (R 4.0.3)
lubridate 1.7.10 2021-02-26 [1] CRAN (R 4.0.3)
magrittr * 2.0.1 2020-11-17 [1] CRAN (R 4.0.2)
MASS 7.3-54 2021-05-03 [1] CRAN (R 4.0.3)
Matrix 1.3-3 2021-05-04 [1] CRAN (R 4.0.3)
mgcv 1.8-35 2021-04-18 [1] CRAN (R 4.0.3)
modelr 0.1.8 2020-05-19 [1] CRAN (R 4.0.2)
multtest 2.46.0 2020-10-27 [1] Bioconductor
munsell 0.5.0 2018-06-12 [1] CRAN (R 4.0.2)
nlme 3.1-152 2021-02-04 [1] CRAN (R 4.0.3)
permute * 0.9-5 2019-03-12 [1] CRAN (R 4.0.2)
phyloseq * 1.34.0 2020-10-27 [1] Bioconductor
pillar 1.6.1 2021-05-16 [1] CRAN (R 4.0.3)
pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.0.2)
plyr 1.8.6 2020-03-03 [1] CRAN (R 4.0.2)
png 0.1-7 2013-12-03 [1] CRAN (R 4.0.2)
prettyunits 1.1.1 2020-01-24 [1] CRAN (R 4.0.2)
progress 1.2.2 2019-05-16 [1] CRAN (R 4.0.2)
purrr * 0.3.4 2020-04-17 [1] CRAN (R 4.0.2)
R6 2.5.0 2020-10-28 [1] CRAN (R 4.0.2)
ragg * 1.1.2 2021-03-17 [1] CRAN (R 4.0.3)
Rcpp 1.0.6 2021-01-15 [1] CRAN (R 4.0.3)
readr * 1.4.0 2020-10-05 [1] CRAN (R 4.0.2)
readxl 1.3.1 2019-03-13 [1] CRAN (R 4.0.2)
reprex 2.0.0 2021-04-02 [1] CRAN (R 4.0.3)
reshape2 1.4.4 2020-04-09 [1] CRAN (R 4.0.2)
rhdf5 2.34.0 2020-10-27 [1] Bioconductor
rhdf5filters 1.2.1 2021-05-03 [1] Bioconductor
Rhdf5lib 1.12.1 2021-01-26 [1] Bioconductor
rlang 0.4.11 2021-04-30 [1] CRAN (R 4.0.3)
rmarkdown * 2.8 2021-05-07 [1] CRAN (R 4.0.3)
rprojroot 2.0.2 2020-11-15 [1] CRAN (R 4.0.2)
rstudioapi 0.13 2020-11-12 [1] CRAN (R 4.0.2)
Rttf2pt1 1.3.8 2020-01-10 [1] CRAN (R 4.0.2)
rvest 1.0.0 2021-03-09 [1] CRAN (R 4.0.3)
S4Vectors * 0.28.1 2020-12-09 [1] Bioconductor
scales * 1.1.1 2020-05-11 [1] CRAN (R 4.0.2)
sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 4.0.2)
stringi 1.6.2 2021-05-17 [1] CRAN (R 4.0.3)
stringr * 1.4.0 2019-02-10 [1] CRAN (R 4.0.2)
survival 3.2-11 2021-04-26 [1] CRAN (R 4.0.3)
svglite * 2.0.0 2021-02-20 [1] CRAN (R 4.0.3)
systemfonts 1.0.2 2021-05-11 [1] CRAN (R 4.0.3)
textshaping 0.3.4 2021-05-11 [1] CRAN (R 4.0.3)
tibble * 3.1.2 2021-05-16 [1] CRAN (R 4.0.3)
tidyr * 1.1.3 2021-03-03 [1] CRAN (R 4.0.3)
tidyselect 1.1.1 2021-04-30 [1] CRAN (R 4.0.3)
tidyverse * 1.3.1 2021-04-15 [1] CRAN (R 4.0.3)
utf8 1.2.1 2021-03-12 [1] CRAN (R 4.0.3)
vctrs 0.3.8 2021-04-29 [1] CRAN (R 4.0.3)
vegan * 2.5-7 2020-11-28 [1] CRAN (R 4.0.3)
viridisLite 0.4.0 2021-04-13 [1] CRAN (R 4.0.3)
webshot 0.5.2 2019-11-22 [1] CRAN (R 4.0.2)
withr 2.4.2 2021-04-18 [1] CRAN (R 4.0.3)
xfun 0.23 2021-05-15 [1] CRAN (R 4.0.3)
xml2 1.3.2 2020-04-23 [1] CRAN (R 4.0.2)
XVector * 0.30.0 2020-10-27 [1] Bioconductor
yaml 2.2.1 2020-02-01 [1] CRAN (R 4.0.2)
zlibbioc 1.36.0 2020-10-27 [1] Bioconductor
[1] /home/angel/R/library
[2] /usr/local/lib/R/site-library
[3] /usr/lib/R/site-library
[4] /usr/lib/R/libraryCallahan BJ, Sankaran K, Fukuyama JA et al. Bioconductor workflow for microbiome data analysis: From raw reads to community analyses. F1000Research 2016;5:1492.